Pre set-up
Import libraries
library(reshape2)
library(ggplot2)
library(plotly)
library(plyr)
library(Rfast)
library(data.table)
library(knitr)
library(rstudioapi)
Import own functions
# Get path current script is in
source_path = dirname(rstudioapi::getSourceEditorContext()$path)
# Source functions required for this script
source(file.path(source_path, 'Beta_pseudo_sim.R', fsep = .Platform$file.sep))
source(file.path(source_path, 'Gaussian_pseudo_sim.R', fsep = .Platform$file.sep))
source(file.path(source_path, 'Create_design.R', fsep = .Platform$file.sep))
source(file.path(source_path, 'Create_miniblock.R', fsep = .Platform$file.sep))
source(file.path(source_path, 'PE_dep_LR_choice_model.R', fsep = .Platform$file.sep))
Parameter set-up
Set up reward space
# Reward lower and upper boundary
reward_space_lb = 1
reward_space_ub = 100
# Set up space for plotting
reward_space_length = 1000
reward_space = seq(reward_space_lb,
reward_space_ub,
length=reward_space_length)
Gaussian parameters and distribution
# Parameters
gaussian_mean = (100 * 1/6) * 3
gaussian_sd = (100 * 1/6) / 3
# Distribution
gaussian_densities = dnorm(reward_space, gaussian_mean, gaussian_sd)
# Normalize density
gaussian_densities = scale(gaussian_densities) - min(scale(gaussian_densities))
Beta parameters and distributions
# Set up "skewdness" parameter for easier change of slope
beta_skewedness = 5
# Parameters for lower end beta and upper end beta
beta_a_le = beta_skewedness/3
beta_b_le = 2*beta_skewedness/3
beta_a_ue = 2*beta_skewedness/3
beta_b_ue = beta_skewedness/3
# Create distributions
beta_densities_le = dbeta(seq(0,1,length=reward_space_length), beta_a_le,beta_b_le)
beta_densities_ue = dbeta(seq(0,1,length=reward_space_length), beta_a_ue,beta_b_ue)
# Normalize densities
beta_densities_le = scale(beta_densities_le) - min(scale(beta_densities_le))
beta_densities_ue = scale(beta_densities_ue) - min(scale(beta_densities_le))
Set rare event threshold
# Set rare events to be 20% of the sampled data
rare_lim = 0.2
Control plots
Plot distributions
# Form data frame from density values
data_densities = data.frame(reward_space,
gaussian_densities,
beta_densities_le,
beta_densities_ue)
# Bring data frame into long format
data_densities = melt(data_densities, id.vars='reward_space')
# Disribution plot
p_dist = ggplot(data=data_densities, aes(x=reward_space, y=value, group=as.factor(variable),
color=as.factor(variable))) +
geom_line(size=1) +
scale_color_brewer(palette='Accent') +
scale_x_continuous(limits = c(1,100)) +
#stat_summary(fun.data = 'mean', geom = 'vline') +
geom_vline(xintercept=c(1/3, 1/2, 2/3)*100, linetype='dashed', alpha=0.5, size=0.1) +
labs(title='Density functions of used distributions',
x='Outcome',
y='Normalized density',
color='Distributions') +
theme(plot.title=element_text(size=15, face='bold', hjust=0.5),
axis.title=element_text(size=12))
p_dist

Test sampling with histograms
# Define parameters
n_samples = 500
data_sampling = data.frame(c(1:n_samples))
colnames(data_sampling) = 'Sample'
# Fill data frame with samples
data_sampling$stim_1 = Beta_pseudo_sim(n_samples, beta_a_le, beta_b_le, 'b_le',
reward_space_lb, reward_space_ub)$outcome
data_sampling$stim_2 = Gaussian_pseudo_sim(n_samples, gaussian_mean, gaussian_sd, 'g_an',
reward_space_lb, reward_space_ub)$outcome
data_sampling$stim_3 = Beta_pseudo_sim(n_samples, beta_a_ue, beta_b_ue, 'b_ue',
reward_space_lb, reward_space_ub)$outcome
data_sampling = melt(data_sampling, id.vars='Sample')
# Plot histograms
p_sampling = ggplot(data=data_sampling, aes(x=value, fill=variable)) +
scale_fill_brewer(palette = 'Accent') +
geom_histogram(binwidth = 1) +
facet_wrap(~variable)
ggplotly(p_sampling)
Simulating data for the choice model
LS0tCnRpdGxlOiAiU2ltdWxhdGlvbiBQRSBkZXBlbmRlbnQgTFIiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCiMgUHJlIHNldC11cAoKIyMjIEltcG9ydCBsaWJyYXJpZXMKYGBge3J9CmxpYnJhcnkocmVzaGFwZTIpCmxpYnJhcnkoZ2dwbG90MikKbGlicmFyeShwbG90bHkpCmxpYnJhcnkocGx5cikKbGlicmFyeShSZmFzdCkKbGlicmFyeShkYXRhLnRhYmxlKQpsaWJyYXJ5KGtuaXRyKQpsaWJyYXJ5KHJzdHVkaW9hcGkpCmBgYAoKCiMjIyBJbXBvcnQgb3duIGZ1bmN0aW9ucwpgYGB7cn0KIyBHZXQgcGF0aCBjdXJyZW50IHNjcmlwdCBpcyBpbgpzb3VyY2VfcGF0aCA9IGRpcm5hbWUocnN0dWRpb2FwaTo6Z2V0U291cmNlRWRpdG9yQ29udGV4dCgpJHBhdGgpCiMgU291cmNlIGZ1bmN0aW9ucyByZXF1aXJlZCBmb3IgdGhpcyBzY3JpcHQKc291cmNlKGZpbGUucGF0aChzb3VyY2VfcGF0aCwgJ0JldGFfcHNldWRvX3NpbS5SJywgZnNlcCA9IC5QbGF0Zm9ybSRmaWxlLnNlcCkpCnNvdXJjZShmaWxlLnBhdGgoc291cmNlX3BhdGgsICdHYXVzc2lhbl9wc2V1ZG9fc2ltLlInLCBmc2VwID0gLlBsYXRmb3JtJGZpbGUuc2VwKSkKc291cmNlKGZpbGUucGF0aChzb3VyY2VfcGF0aCwgJ0NyZWF0ZV9kZXNpZ24uUicsIGZzZXAgPSAuUGxhdGZvcm0kZmlsZS5zZXApKQpzb3VyY2UoZmlsZS5wYXRoKHNvdXJjZV9wYXRoLCAnQ3JlYXRlX21pbmlibG9jay5SJywgZnNlcCA9IC5QbGF0Zm9ybSRmaWxlLnNlcCkpCnNvdXJjZShmaWxlLnBhdGgoc291cmNlX3BhdGgsICdQRV9kZXBfTFJfY2hvaWNlX21vZGVsLlInLCBmc2VwID0gLlBsYXRmb3JtJGZpbGUuc2VwKSkKYGBgCgojIFBhcmFtZXRlciBzZXQtdXAKCiMjIyBTZXQgdXAgcmV3YXJkIHNwYWNlCmBgYHtyfQojIFJld2FyZCBsb3dlciBhbmQgdXBwZXIgYm91bmRhcnkKcmV3YXJkX3NwYWNlX2xiID0gMQpyZXdhcmRfc3BhY2VfdWIgPSAxMDAKCiMgU2V0IHVwIHNwYWNlIGZvciBwbG90dGluZwpyZXdhcmRfc3BhY2VfbGVuZ3RoID0gMTAwMApyZXdhcmRfc3BhY2UgPSBzZXEocmV3YXJkX3NwYWNlX2xiLCAKICAgICAgICAgICAgICAgICAgIHJld2FyZF9zcGFjZV91YiwKICAgICAgICAgICAgICAgICAgIGxlbmd0aD1yZXdhcmRfc3BhY2VfbGVuZ3RoKQpgYGAKCiMjIyBHYXVzc2lhbiBwYXJhbWV0ZXJzIGFuZCBkaXN0cmlidXRpb24KYGBge3J9CiMgUGFyYW1ldGVycwpnYXVzc2lhbl9tZWFuID0gKDEwMCAqIDEvNikgKiAzCmdhdXNzaWFuX3NkID0gKDEwMCAqIDEvNikgLyAzCgojIERpc3RyaWJ1dGlvbgpnYXVzc2lhbl9kZW5zaXRpZXMgPSBkbm9ybShyZXdhcmRfc3BhY2UsIGdhdXNzaWFuX21lYW4sIGdhdXNzaWFuX3NkKQojIE5vcm1hbGl6ZSBkZW5zaXR5CmdhdXNzaWFuX2RlbnNpdGllcyA9IHNjYWxlKGdhdXNzaWFuX2RlbnNpdGllcykgLSBtaW4oc2NhbGUoZ2F1c3NpYW5fZGVuc2l0aWVzKSkKYGBgCgojIyMgQmV0YSBwYXJhbWV0ZXJzIGFuZCBkaXN0cmlidXRpb25zCmBgYHtyfQojIFNldCB1cCAic2tld2RuZXNzIiBwYXJhbWV0ZXIgZm9yIGVhc2llciBjaGFuZ2Ugb2Ygc2xvcGUKYmV0YV9za2V3ZWRuZXNzID0gNQojIFBhcmFtZXRlcnMgZm9yIGxvd2VyIGVuZCBiZXRhIGFuZCB1cHBlciBlbmQgYmV0YQpiZXRhX2FfbGUgPSBiZXRhX3NrZXdlZG5lc3MvMwpiZXRhX2JfbGUgPSAyKmJldGFfc2tld2VkbmVzcy8zCmJldGFfYV91ZSA9IDIqYmV0YV9za2V3ZWRuZXNzLzMKYmV0YV9iX3VlID0gYmV0YV9za2V3ZWRuZXNzLzMKIyBDcmVhdGUgZGlzdHJpYnV0aW9ucwpiZXRhX2RlbnNpdGllc19sZSA9IGRiZXRhKHNlcSgwLDEsbGVuZ3RoPXJld2FyZF9zcGFjZV9sZW5ndGgpLCBiZXRhX2FfbGUsYmV0YV9iX2xlKQpiZXRhX2RlbnNpdGllc191ZSA9IGRiZXRhKHNlcSgwLDEsbGVuZ3RoPXJld2FyZF9zcGFjZV9sZW5ndGgpLCBiZXRhX2FfdWUsYmV0YV9iX3VlKQojIE5vcm1hbGl6ZSBkZW5zaXRpZXMKYmV0YV9kZW5zaXRpZXNfbGUgPSBzY2FsZShiZXRhX2RlbnNpdGllc19sZSkgLSBtaW4oc2NhbGUoYmV0YV9kZW5zaXRpZXNfbGUpKQpiZXRhX2RlbnNpdGllc191ZSA9IHNjYWxlKGJldGFfZGVuc2l0aWVzX3VlKSAtIG1pbihzY2FsZShiZXRhX2RlbnNpdGllc19sZSkpCmBgYAoKIyMjIFNldCByYXJlIGV2ZW50IHRocmVzaG9sZApgYGB7cn0KIyBTZXQgcmFyZSBldmVudHMgdG8gYmUgMjAlIG9mIHRoZSBzYW1wbGVkIGRhdGEKcmFyZV9saW0gPSAwLjIKYGBgCgotLS0KCiMgQ29udHJvbCBwbG90cwoKIyMjIFBsb3QgZGlzdHJpYnV0aW9ucwpgYGB7ciBvdXQud2lkdGg9JzEwMCUnLCB3YXJuaW5nPUZBTFNFfQojIEZvcm0gZGF0YSBmcmFtZSBmcm9tIGRlbnNpdHkgdmFsdWVzCmRhdGFfZGVuc2l0aWVzID0gZGF0YS5mcmFtZShyZXdhcmRfc3BhY2UsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgZ2F1c3NpYW5fZGVuc2l0aWVzLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgYmV0YV9kZW5zaXRpZXNfbGUsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICBiZXRhX2RlbnNpdGllc191ZSkKCiMgQnJpbmcgZGF0YSBmcmFtZSBpbnRvIGxvbmcgZm9ybWF0CmRhdGFfZGVuc2l0aWVzID0gbWVsdChkYXRhX2RlbnNpdGllcywgaWQudmFycz0ncmV3YXJkX3NwYWNlJykKCiMgRGlzcmlidXRpb24gcGxvdApwX2Rpc3QgPSBnZ3Bsb3QoZGF0YT1kYXRhX2RlbnNpdGllcywgYWVzKHg9cmV3YXJkX3NwYWNlLCB5PXZhbHVlLCBncm91cD1hcy5mYWN0b3IodmFyaWFibGUpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBjb2xvcj1hcy5mYWN0b3IodmFyaWFibGUpKSkgKwogIGdlb21fbGluZShzaXplPTEpICsKICBzY2FsZV9jb2xvcl9icmV3ZXIocGFsZXR0ZT0nQWNjZW50JykgKwogIHNjYWxlX3hfY29udGludW91cyhsaW1pdHMgPSBjKDEsMTAwKSkgKwogICNzdGF0X3N1bW1hcnkoZnVuLmRhdGEgPSAnbWVhbicsIGdlb20gPSAndmxpbmUnKSArCiAgZ2VvbV92bGluZSh4aW50ZXJjZXB0PWMoMS8zLCAxLzIsIDIvMykqMTAwLCBsaW5ldHlwZT0nZGFzaGVkJywgYWxwaGE9MC41LCBzaXplPTAuMSkgKwogIGxhYnModGl0bGU9J0RlbnNpdHkgZnVuY3Rpb25zIG9mIHVzZWQgZGlzdHJpYnV0aW9ucycsCiAgICAgICB4PSdPdXRjb21lJywKICAgICAgIHk9J05vcm1hbGl6ZWQgZGVuc2l0eScsCiAgICAgICBjb2xvcj0nRGlzdHJpYnV0aW9ucycpICsKICB0aGVtZShwbG90LnRpdGxlPWVsZW1lbnRfdGV4dChzaXplPTE1LCBmYWNlPSdib2xkJywgaGp1c3Q9MC41KSwKICAgICAgICBheGlzLnRpdGxlPWVsZW1lbnRfdGV4dChzaXplPTEyKSkKCnBfZGlzdApgYGAKCiMjIyBUZXN0IHNhbXBsaW5nIHdpdGggaGlzdG9ncmFtcwpgYGB7ciBvdXQud2lkdGg9JzEwMCUnLCB3YXJuaW5nPUZBTFNFfQojIERlZmluZSBwYXJhbWV0ZXJzCm5fc2FtcGxlcyA9IDUwMApkYXRhX3NhbXBsaW5nID0gZGF0YS5mcmFtZShjKDE6bl9zYW1wbGVzKSkKY29sbmFtZXMoZGF0YV9zYW1wbGluZykgPSAnU2FtcGxlJwoKIyBGaWxsIGRhdGEgZnJhbWUgd2l0aCBzYW1wbGVzCmRhdGFfc2FtcGxpbmckc3RpbV8xID0gQmV0YV9wc2V1ZG9fc2ltKG5fc2FtcGxlcywgYmV0YV9hX2xlLCBiZXRhX2JfbGUsICdiX2xlJywKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgcmV3YXJkX3NwYWNlX2xiLCByZXdhcmRfc3BhY2VfdWIpJG91dGNvbWUKZGF0YV9zYW1wbGluZyRzdGltXzIgPSBHYXVzc2lhbl9wc2V1ZG9fc2ltKG5fc2FtcGxlcywgZ2F1c3NpYW5fbWVhbiwgZ2F1c3NpYW5fc2QsICdnX2FuJywKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHJld2FyZF9zcGFjZV9sYiwgcmV3YXJkX3NwYWNlX3ViKSRvdXRjb21lCmRhdGFfc2FtcGxpbmckc3RpbV8zID0gQmV0YV9wc2V1ZG9fc2ltKG5fc2FtcGxlcywgYmV0YV9hX3VlLCBiZXRhX2JfdWUsICdiX3VlJywKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgcmV3YXJkX3NwYWNlX2xiLCByZXdhcmRfc3BhY2VfdWIpJG91dGNvbWUKZGF0YV9zYW1wbGluZyA9IG1lbHQoZGF0YV9zYW1wbGluZywgaWQudmFycz0nU2FtcGxlJykKCiMgUGxvdCBoaXN0b2dyYW1zCnBfc2FtcGxpbmcgPSBnZ3Bsb3QoZGF0YT1kYXRhX3NhbXBsaW5nLCBhZXMoeD12YWx1ZSwgZmlsbD12YXJpYWJsZSkpICsKICBzY2FsZV9maWxsX2JyZXdlcihwYWxldHRlID0gJ0FjY2VudCcpICsKICBnZW9tX2hpc3RvZ3JhbShiaW53aWR0aCA9IDEpICsKICBmYWNldF93cmFwKH52YXJpYWJsZSkKCmdncGxvdGx5KHBfc2FtcGxpbmcpCmBgYAoKLS0tCgojIyMgU2ltdWxhdGluZyBkYXRhIGZvciB0aGUgY2hvaWNlIG1vZGVsCgpgYGB7ciBlY2hvPUZBTFNFfQojIERlZmluZSBzaW11bGF0aW9uIHBhcmFtZXRlcnMKbl9zdWJqZWN0cyA9IDQwCnN1YmplY3RzID0gYygxOm5fc3ViamVjdHMpCm5fbWluaWJsb2NrcyA9IDI1CmFscGhhMCA9IDAuMwphbHBoYTEgPSAwLjYKdGVtcGVyYXR1cmUgPSA3CgojIFNldCB1cCBtYXRyaXggdG8gc3RvcmUgZGF0YSBpbiAodG8gYXBwZW5kIHN1YmplY3RzKQpkYXRhX21vZGVsID0gbWF0cml4KDAsMCwxNykKCiMgTG9vcCBvdmVyIGVhY2ggc3ViamVjdApmb3Ioc3ViX2NvdW50IGluIHN1YmplY3RzKXsKICAjIFNpbXVsYXRlIGRhdGEgZm9yIHN1YmplY3QKICBkZXNpZ24gPSBDcmVhdGVfZGVzaWduKG5fbWluaWJsb2NrcywKICAgICAgICAgICAgICAgICAgICAgICAgIGJldGFfYV9sZSwKICAgICAgICAgICAgICAgICAgICAgICAgIGJldGFfYl9sZSwKICAgICAgICAgICAgICAgICAgICAgICAgIGdhdXNzaWFuX21lYW4sCiAgICAgICAgICAgICAgICAgICAgICAgICBnYXVzc2lhbl9zZCwKICAgICAgICAgICAgICAgICAgICAgICAgIGJldGFfYV91ZSwKICAgICAgICAgICAgICAgICAgICAgICAgIGJldGFfYl91ZSwKICAgICAgICAgICAgICAgICAgICAgICAgIHJld2FyZF9zcGFjZV9sYiwKICAgICAgICAgICAgICAgICAgICAgICAgIHJld2FyZF9zcGFjZV91YikKICAKICAjIENyZWF0ZSBtYXRyaXggdG8gc3RvcmUgZGF0YSBmb3IgZWFjaCBzdWJqZWN0CiAgcmVzdWx0cyA9IG1hdHJpeCgwLG5yb3coZGVzaWduKSxuY29sKGRhdGFfbW9kZWwpKQogIHJlc3VsdHNbLDFdID0gc3ViamVjdHNbc3ViX2NvdW50XQogIHJlc3VsdHNbLDJdID0gYygxOm5yb3coZGVzaWduKSkKICByZXN1bHRzWywzXSA9IGRlc2lnbiRvcHRpb25fbGVmdAogIHJlc3VsdHNbLDRdID0gZGVzaWduJG9wdGlvbl9yaWdodAogIHJlc3VsdHNbLDVdID0gZGVzaWduJHJld2FyZF9zdGltXzEKICByZXN1bHRzWyw2XSA9IGRlc2lnbiRyZXdhcmRfc3RpbV8yCiAgcmVzdWx0c1ssN10gPSBkZXNpZ24kY29tcF9udW1iZXIKCgogICMgUnVuIG1vZGVsIG92ZXIgc2ltdWxhdGVkIGRhdGEKICBzaW0gPSBQRV9kZXBfTFJfY2hvaWNlX21vZGVsKGRlc2lnbiwgYWxwaGEwLCBhbHBoYTEsIHRlbXBlcmF0dXJlKQogICMgU2F2ZSBtb2RlbCB2YWx1ZXMsIFBFIGFuZCBmUEUgZm9yIGVjaCB0cmlhbCBhbmQgc3ViamVjdCBpbnRvIG1hdHJpeAogIHJlc3VsdHNbLDhdID0gc2ltJGNob2ljZSRjaG9pY2UKICByZXN1bHRzWyw5OjExXSA9IGFzLm1hdHJpeChzaW0kdmFsdWVzKQogIHJlc3VsdHNbLDEyOjE0XSA9IGFzLm1hdHJpeChzaW0kUEUpCiAgcmVzdWx0c1ssMTU6MTddID0gYXMubWF0cml4KHNpbSRmUEUpCgogICMgQXBwZW5kIGFsbCBzdWJqZWN0cwogIGRhdGFfbW9kZWwgPSByYmluZChkYXRhX21vZGVsLCByZXN1bHRzKQp9CgojIENvbnZlcnQgdG8gZGF0YSBmcmFtZSAoSG93IG91ciBkYXRhIHdpbGwgbG9vayBsaWtlIGluIHRoZSBzdHVkeSkKZGF0YV9tb2RlbCA9IGRhdGEuZnJhbWUoZGF0YV9tb2RlbCkKY29sbmFtZXMoZGF0YV9tb2RlbCkgPSBjKCdzdWJfaWQnLCAndHJpYWwnLCAnc3RpbV8xJywgJ3N0aW1fMicsCiAgICAgICAgICAgICAgICAgICAgICAgICAncmV3YXJkX3N0aW1fMScsICdyZXdhcmRfc3RpbV8yJywgJ2NvbXBfbnVtYmVyJywgJ2Nob2ljZScsCiAgICAgICAgICAgICAgICAgICAgICAgICd2X3N0aW1fMScsICd2X3N0aW1fMicsICd2X3N0aW1fMycsCiAgICAgICAgICAgICAgICAgICAgICAgICdwZV9zdGltXzEnLCAncGVfc3RpbV8yJywgJ3BlX3N0aW1fMycsCiAgICAgICAgICAgICAgICAgICAgICAgICdmcGVfc3RpbV8xJywgJ2ZwZV9zdGltXzInLCAnZnBlX3N0aW1fMycpCiMgU29ydCBhZnRlciBwYXJ0aWNpcGFudHMKZGF0YV9tb2RlbCA9IGRhdGFfbW9kZWxbb3JkZXIoZGF0YV9tb2RlbCRzdWJfaWQpLF0KYGBgCgo=